## Section 3.2.2 of the pre-analysis plan: "Overall Precision: Expected
## Distance Between Ranks Drawn from a Given Posterior Distribution"
##
## Kevin Quinn
## University of Michigan
##
## 8/7/2019
##

library(TopKLists)
library(coda)

set.seed(44919)

sink(file="MM3.2.2.output.txt")



########################################################################
## START a1-b1 comparisons
########################################################################


## load MCMC output
load("MM3.2.a1.M1.out.Rda")
a1.M1.out <- a1.M1.out[sample(1:nrow(a1.M1.out), nrow(a1.M1.out),
                              replace=FALSE), ]
load("MM3.2.a1.M2.out.Rda")
a1.M2.out <- a1.M2.out[sample(1:nrow(a1.M2.out), nrow(a1.M2.out),
                              replace=FALSE), ]
load("MM3.2.a1.pairwise.out.Rda")
a1.pairwise.out <- a1.pairwise.out[sample(1:nrow(a1.pairwise.out),
                                          nrow(a1.pairwise.out),
                              replace=FALSE), ]

load("MM3.2.b1.M1.out.Rda")
b1.M1.out <- b1.M1.out[sample(1:nrow(b1.M1.out), nrow(b1.M1.out),
                              replace=FALSE), ]

load("MM3.2.b1.M2.out.Rda")
b1.M2.out <- b1.M2.out[sample(1:nrow(b1.M2.out), nrow(b1.M2.out),
                              replace=FALSE), ]
load("MM3.2.b1.pairwise.out.Rda")
b1.pairwise.out <- b1.pairwise.out[sample(1:nrow(b1.pairwise.out),
                                          nrow(b1.pairwise.out),
                              replace=FALSE), ]


## make MCMC output comparable
a1.M1.out <- a1.M1.out[,1:48]
a1.M2.out <- a1.M2.out[,1:48]
b1.M1.out <- b1.M1.out[,1:48]
b1.M2.out <- b1.M2.out[,1:48]

colnames(a1.M1.out) <- gsub("as.factor\\(photoID\\)", "", colnames(a1.M1.out))
colnames(a1.M2.out) <- gsub("as.factor\\(photoID\\)", "", colnames(a1.M2.out))
colnames(b1.M1.out) <- gsub("as.factor\\(photoID\\)", "", colnames(b1.M1.out))
colnames(b1.M2.out) <- gsub("as.factor\\(photoID\\)", "", colnames(b1.M2.out))
colnames(a1.pairwise.out) <- gsub("theta\\.", "", colnames(a1.pairwise.out))
colnames(b1.pairwise.out) <- gsub("theta\\.", "", colnames(b1.pairwise.out))

for (j in 2:48){
    a1.M1.out[,j] <- a1.M1.out[,1] + a1.M1.out[,j]
    a1.M2.out[,j] <- a1.M2.out[,1] + a1.M2.out[,j]
    b1.M1.out[,j] <- b1.M1.out[,1] + b1.M1.out[,j]
    b1.M2.out[,j] <- b1.M2.out[,1] + b1.M2.out[,j]
}
colnames(a1.M1.out)[1] <- colnames(a1.pairwise.out)[1]
colnames(a1.M2.out)[1] <- colnames(a1.pairwise.out)[1]
colnames(b1.M1.out)[1] <- colnames(a1.pairwise.out)[1]
colnames(b1.M2.out)[1] <- colnames(a1.pairwise.out)[1]


M <- nrow(a1.M1.out)

## permute the MCMC output
perm.inds <- sample(1:M, size=M, replace=FALSE)
a1.M1.out.perm <- a1.M1.out[perm.inds,]
perm.inds <- sample(1:M, size=M, replace=FALSE)
a1.M2.out.perm <- a1.M2.out[perm.inds,]
perm.inds <- sample(1:M, size=M, replace=FALSE)
b1.M1.out.perm <- b1.M1.out[perm.inds,]
perm.inds <- sample(1:M, size=M, replace=FALSE)
b1.M2.out.perm <- b1.M2.out[perm.inds,]
perm.inds <- sample(1:M, size=M, replace=FALSE)
a1.pairwise.out.perm <- a1.pairwise.out[perm.inds,]
perm.inds <- sample(1:M, size=M, replace=FALSE)
b1.pairwise.out.perm <- b1.pairwise.out[perm.inds,]


## storage objects
D.spear.a1.M1.holder <- rep(NA, M)
D.spear.a1.M2.holder <- rep(NA, M)
D.spear.b1.M1.holder <- rep(NA, M)
D.spear.b1.M2.holder <- rep(NA, M)
D.spear.a1.pair.holder <- rep(NA, M)
D.spear.b1.pair.holder <- rep(NA, M)

D.kend.a1.M1.holder <- rep(NA, M)
D.kend.a1.M2.holder <- rep(NA, M)
D.kend.b1.M1.holder <- rep(NA, M)
D.kend.b1.M2.holder <- rep(NA, M)
D.kend.a1.pair.holder <- rep(NA, M)
D.kend.b1.pair.holder <- rep(NA, M)


## M1 vs. pairwise
for (iter in 1:M){


    D.spear.a1.M1.holder[iter] <- Spearman(rank(a1.M1.out[iter,]),
                                           rank(a1.M1.out.perm[iter,]),
                                           48, 48)
    D.spear.a1.M2.holder[iter] <- Spearman(rank(a1.M2.out[iter,]),
                                           rank(a1.M2.out.perm[iter,]),
                                           48, 48)
    D.spear.b1.M1.holder[iter] <- Spearman(rank(b1.M1.out[iter,]),
                                           rank(b1.M1.out.perm[iter,]),
                                           48, 48)
    D.spear.b1.M2.holder[iter] <- Spearman(rank(b1.M2.out[iter,]),
                                           rank(b1.M2.out.perm[iter,]),
                                           48, 48)
    D.spear.a1.pair.holder[iter] <- Spearman(rank(a1.pairwise.out[iter,]),
                                             rank(a1.pairwise.out.perm[iter,]),
                                             48, 48)
    D.spear.b1.pair.holder[iter] <- Spearman(rank(b1.pairwise.out[iter,]),
                                             rank(b1.pairwise.out.perm[iter,]),
                                             48, 48)

    
    D.kend.a1.M1.holder[iter] <- Kendall2Lists(rank(a1.M1.out[iter,]),
                                               rank(a1.M1.out.perm[iter,]),
                                               48, 48, 48)
    D.kend.a1.M2.holder[iter] <- Kendall2Lists(rank(a1.M2.out[iter,]),
                                               rank(a1.M2.out.perm[iter,]),
                                               48, 48, 48)
    D.kend.b1.M1.holder[iter] <- Kendall2Lists(rank(b1.M1.out[iter,]),
                                               rank(b1.M1.out.perm[iter,]),
                                               48, 48, 48)
    D.kend.b1.M2.holder[iter] <- Kendall2Lists(rank(b1.M2.out[iter,]),
                                               rank(b1.M2.out.perm[iter,]),
                                               48, 48, 48)
    D.kend.a1.pair.holder[iter] <- Kendall2Lists(rank(a1.pairwise.out[iter,]),
                                                 rank(a1.pairwise.out.perm[iter,]),
                                                 48, 48, 48)
    D.kend.b1.pair.holder[iter] <- Kendall2Lists(rank(b1.pairwise.out[iter,]),
                                                 rank(b1.pairwise.out.perm[iter,]),
                                                 48, 48, 48)
    
}



cat("\n\n@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@\n")
cat("a1 - b1 precision (coder race correlated with photos)\n")
cat("-----------------------------------------------------------\n")
cat("                        D.spear                 D.kend\n")
cat("-----------------------------------------------------------\n")
cat("a1.M1                   ", round(mean(D.spear.a1.M1.holder), 3),
    "               ",
    round(mean(D.kend.a1.M1.holder), 3), "\n")
cat("b1.M1                   ", round(mean(D.spear.b1.M1.holder), 3),
    "               ",
    round(mean(D.kend.b1.M1.holder), 3), "\n")
cat("a1.M2                   ", round(mean(D.spear.a1.M2.holder), 3),
    "               ",
    round(mean(D.kend.a1.M2.holder), 3), "\n")
cat("b1.M2                   ", round(mean(D.spear.b1.M2.holder), 3),
    "               ",
    round(mean(D.kend.b1.M2.holder), 3), "\n")
cat("a1.pairwise             ", round(mean(D.spear.a1.pair.holder), 3),
    "               ",
    round(mean(D.kend.a1.pair.holder), 3), "\n")
cat("b1.pairwise             ", round(mean(D.spear.b1.pair.holder), 3),
    "               ",
    round(mean(D.kend.b1.pair.holder), 3), "\n")

cat("@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@\n\n")


D.hat.data <- data.frame(D.spear.a1.M1=D.spear.a1.M1.holder,
                         D.spear.a1.M2=D.spear.a1.M2.holder,
                         D.spear.b1.M1=D.spear.b1.M1.holder,
                         D.spear.b1.M2=D.spear.b1.M2.holder,
                         D.spear.a1.pair=D.spear.a1.pair.holder,
                         D.spear.b1.pair=D.spear.b1.pair.holder,
                         D.kend.a1.M1=D.kend.a1.M1.holder,
                         D.kend.a1.M2=D.kend.a1.M2.holder,
                         D.kend.b1.M1=D.kend.b1.M1.holder,
                         D.kend.b1.M2=D.kend.b1.M2.holder,
                         D.kend.a1.pair=D.kend.a1.pair.holder,
                         D.kend.b1.pair=D.kend.b1.pair.holder)


save(D.hat.data, file="MM3.2.2.D.hat.a1-b1.data.Rda")
                         
                         



########################################################################
## END a1-b1 comparisons
########################################################################











########################################################################
## START a2-b2 comparisons
########################################################################


## load MCMC output
load("MM3.2.a2.M1.out.Rda")
a2.M1.out <- a2.M1.out[sample(1:nrow(a2.M1.out), nrow(a2.M1.out),
                              replace=FALSE), ]
load("MM3.2.a2.M2.out.Rda")
a2.M2.out <- a2.M2.out[sample(1:nrow(a2.M2.out), nrow(a2.M2.out),
                              replace=FALSE), ]
load("MM3.2.a2.pairwise.out.Rda")
a2.pairwise.out <- a2.pairwise.out[sample(1:nrow(a2.pairwise.out),
                                          nrow(a2.pairwise.out),
                              replace=FALSE), ]

load("MM3.2.b2.M1.out.Rda")
b2.M1.out <- b2.M1.out[sample(1:nrow(b2.M1.out), nrow(b2.M1.out),
                              replace=FALSE), ]

load("MM3.2.b2.M2.out.Rda")
b2.M2.out <- b2.M2.out[sample(1:nrow(b2.M2.out), nrow(b2.M2.out),
                              replace=FALSE), ]
load("MM3.2.b2.pairwise.out.Rda")
b2.pairwise.out <- b2.pairwise.out[sample(1:nrow(b2.pairwise.out),
                                          nrow(b2.pairwise.out),
                              replace=FALSE), ]


## make MCMC output comparable
a2.M1.out <- a2.M1.out[,1:48]
a2.M2.out <- a2.M2.out[,1:48]
b2.M1.out <- b2.M1.out[,1:48]
b2.M2.out <- b2.M2.out[,1:48]

colnames(a2.M1.out) <- gsub("as.factor\\(photoID\\)", "", colnames(a2.M1.out))
colnames(a2.M2.out) <- gsub("as.factor\\(photoID\\)", "", colnames(a2.M2.out))
colnames(b2.M1.out) <- gsub("as.factor\\(photoID\\)", "", colnames(b2.M1.out))
colnames(b2.M2.out) <- gsub("as.factor\\(photoID\\)", "", colnames(b2.M2.out))
colnames(a2.pairwise.out) <- gsub("theta\\.", "", colnames(a2.pairwise.out))
colnames(b2.pairwise.out) <- gsub("theta\\.", "", colnames(b2.pairwise.out))

for (j in 2:48){
    a2.M1.out[,j] <- a2.M1.out[,1] + a2.M1.out[,j]
    a2.M2.out[,j] <- a2.M2.out[,1] + a2.M2.out[,j]
    b2.M1.out[,j] <- b2.M1.out[,1] + b2.M1.out[,j]
    b2.M2.out[,j] <- b2.M2.out[,1] + b2.M2.out[,j]
}
colnames(a2.M1.out)[1] <- colnames(a2.pairwise.out)[1]
colnames(a2.M2.out)[1] <- colnames(a2.pairwise.out)[1]
colnames(b2.M1.out)[1] <- colnames(a2.pairwise.out)[1]
colnames(b2.M2.out)[1] <- colnames(a2.pairwise.out)[1]


M <- nrow(a2.M1.out)

## permute the MCMC output
perm.inds <- sample(1:M, size=M, replace=FALSE)
a2.M1.out.perm <- a2.M1.out[perm.inds,]
perm.inds <- sample(1:M, size=M, replace=FALSE)
a2.M2.out.perm <- a2.M2.out[perm.inds,]
perm.inds <- sample(1:M, size=M, replace=FALSE)
b2.M1.out.perm <- b2.M1.out[perm.inds,]
perm.inds <- sample(1:M, size=M, replace=FALSE)
b2.M2.out.perm <- b2.M2.out[perm.inds,]
perm.inds <- sample(1:M, size=M, replace=FALSE)
a2.pairwise.out.perm <- a2.pairwise.out[perm.inds,]
perm.inds <- sample(1:M, size=M, replace=FALSE)
b2.pairwise.out.perm <- b2.pairwise.out[perm.inds,]


## storage objects
D.spear.a2.M1.holder <- rep(NA, M)
D.spear.a2.M2.holder <- rep(NA, M)
D.spear.b2.M1.holder <- rep(NA, M)
D.spear.b2.M2.holder <- rep(NA, M)
D.spear.a2.pair.holder <- rep(NA, M)
D.spear.b2.pair.holder <- rep(NA, M)

D.kend.a2.M1.holder <- rep(NA, M)
D.kend.a2.M2.holder <- rep(NA, M)
D.kend.b2.M1.holder <- rep(NA, M)
D.kend.b2.M2.holder <- rep(NA, M)
D.kend.a2.pair.holder <- rep(NA, M)
D.kend.b2.pair.holder <- rep(NA, M)


## M1 vs. pairwise
for (iter in 1:M){


    D.spear.a2.M1.holder[iter] <- Spearman(rank(a2.M1.out[iter,]),
                                           rank(a2.M1.out.perm[iter,]),
                                           48, 48)
    D.spear.a2.M2.holder[iter] <- Spearman(rank(a2.M2.out[iter,]),
                                           rank(a2.M2.out.perm[iter,]),
                                           48, 48)
    D.spear.b2.M1.holder[iter] <- Spearman(rank(b2.M1.out[iter,]),
                                           rank(b2.M1.out.perm[iter,]),
                                           48, 48)
    D.spear.b2.M2.holder[iter] <- Spearman(rank(b2.M2.out[iter,]),
                                           rank(b2.M2.out.perm[iter,]),
                                           48, 48)
    D.spear.a2.pair.holder[iter] <- Spearman(rank(a2.pairwise.out[iter,]),
                                             rank(a2.pairwise.out.perm[iter,]),
                                             48, 48)
    D.spear.b2.pair.holder[iter] <- Spearman(rank(b2.pairwise.out[iter,]),
                                             rank(b2.pairwise.out.perm[iter,]),
                                             48, 48)

    
    D.kend.a2.M1.holder[iter] <- Kendall2Lists(rank(a2.M1.out[iter,]),
                                               rank(a2.M1.out.perm[iter,]),
                                               48, 48, 48)
    D.kend.a2.M2.holder[iter] <- Kendall2Lists(rank(a2.M2.out[iter,]),
                                               rank(a2.M2.out.perm[iter,]),
                                               48, 48, 48)
    D.kend.b2.M1.holder[iter] <- Kendall2Lists(rank(b2.M1.out[iter,]),
                                               rank(b2.M1.out.perm[iter,]),
                                               48, 48, 48)
    D.kend.b2.M2.holder[iter] <- Kendall2Lists(rank(b2.M2.out[iter,]),
                                               rank(b2.M2.out.perm[iter,]),
                                               48, 48, 48)
    D.kend.a2.pair.holder[iter] <- Kendall2Lists(rank(a2.pairwise.out[iter,]),
                                                 rank(a2.pairwise.out.perm[iter,]),
                                                 48, 48, 48)
    D.kend.b2.pair.holder[iter] <- Kendall2Lists(rank(b2.pairwise.out[iter,]),
                                                 rank(b2.pairwise.out.perm[iter,]),
                                                 48, 48, 48)
    
}



cat("\n\n@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@\n")
cat("a2 - b2 precision (distractor race correlated with photos)\n")
cat("-----------------------------------------------------------\n")
cat("                        D.spear                 D.kend\n")
cat("-----------------------------------------------------------\n")
cat("a2.M1                   ", round(mean(D.spear.a2.M1.holder), 3),
    "               ",
    round(mean(D.kend.a2.M1.holder), 3), "\n")
cat("b2.M1                   ", round(mean(D.spear.b2.M1.holder), 3),
    "               ",
    round(mean(D.kend.b2.M1.holder), 3), "\n")
cat("a2.M2                   ", round(mean(D.spear.a2.M2.holder), 3),
    "               ",
    round(mean(D.kend.a2.M2.holder), 3), "\n")
cat("b2.M2                   ", round(mean(D.spear.b2.M2.holder), 3),
    "               ",
    round(mean(D.kend.b2.M2.holder), 3), "\n")
cat("a2.pairwise             ", round(mean(D.spear.a2.pair.holder), 3),
    "               ",
    round(mean(D.kend.a2.pair.holder), 3), "\n")
cat("b2.pairwise             ", round(mean(D.spear.b2.pair.holder), 3),
    "               ",
    round(mean(D.kend.b2.pair.holder), 3), "\n")

cat("@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@\n\n")


D.hat.data <- data.frame(D.spear.a2.M1=D.spear.a2.M1.holder,
                         D.spear.a2.M2=D.spear.a2.M2.holder,
                         D.spear.b2.M1=D.spear.b2.M1.holder,
                         D.spear.b2.M2=D.spear.b2.M2.holder,
                         D.spear.a2.pair=D.spear.a2.pair.holder,
                         D.spear.b2.pair=D.spear.b2.pair.holder,
                         D.kend.a2.M1=D.kend.a2.M1.holder,
                         D.kend.a2.M2=D.kend.a2.M2.holder,
                         D.kend.b2.M1=D.kend.b2.M1.holder,
                         D.kend.b2.M2=D.kend.b2.M2.holder,
                         D.kend.a2.pair=D.kend.a2.pair.holder,
                         D.kend.b2.pair=D.kend.b2.pair.holder)


save(D.hat.data, file="MM3.2.2.D.hat.a2-b2.data.Rda")
                         
                         



########################################################################
## END a2-b2 comparisons
########################################################################











########################################################################
## START a3-b3 comparisons
########################################################################


## load MCMC output
load("MM3.2.a3.M1.out.Rda")
a3.M1.out <- a3.M1.out[sample(1:nrow(a3.M1.out), nrow(a3.M1.out),
                              replace=FALSE), ]
load("MM3.2.a3.M2.out.Rda")
a3.M2.out <- a3.M2.out[sample(1:nrow(a3.M2.out), nrow(a3.M2.out),
                              replace=FALSE), ]
load("MM3.2.a3.pairwise.out.Rda")
a3.pairwise.out <- a3.pairwise.out[sample(1:nrow(a3.pairwise.out),
                                          nrow(a3.pairwise.out),
                              replace=FALSE), ]

load("MM3.2.b3.M1.out.Rda")
b3.M1.out <- b3.M1.out[sample(1:nrow(b3.M1.out), nrow(b3.M1.out),
                              replace=FALSE), ]

load("MM3.2.b3.M2.out.Rda")
b3.M2.out <- b3.M2.out[sample(1:nrow(b3.M2.out), nrow(b3.M2.out),
                              replace=FALSE), ]
load("MM3.2.b3.pairwise.out.Rda")
b3.pairwise.out <- b3.pairwise.out[sample(1:nrow(b3.pairwise.out),
                                          nrow(b3.pairwise.out),
                              replace=FALSE), ]


## make MCMC output comparable
a3.M1.out <- a3.M1.out[,1:48]
a3.M2.out <- a3.M2.out[,1:48]
b3.M1.out <- b3.M1.out[,1:48]
b3.M2.out <- b3.M2.out[,1:48]

colnames(a3.M1.out) <- gsub("as.factor\\(photoID\\)", "", colnames(a3.M1.out))
colnames(a3.M2.out) <- gsub("as.factor\\(photoID\\)", "", colnames(a3.M2.out))
colnames(b3.M1.out) <- gsub("as.factor\\(photoID\\)", "", colnames(b3.M1.out))
colnames(b3.M2.out) <- gsub("as.factor\\(photoID\\)", "", colnames(b3.M2.out))
colnames(a3.pairwise.out) <- gsub("theta\\.", "", colnames(a3.pairwise.out))
colnames(b3.pairwise.out) <- gsub("theta\\.", "", colnames(b3.pairwise.out))

for (j in 2:48){
    a3.M1.out[,j] <- a3.M1.out[,1] + a3.M1.out[,j]
    a3.M2.out[,j] <- a3.M2.out[,1] + a3.M2.out[,j]
    b3.M1.out[,j] <- b3.M1.out[,1] + b3.M1.out[,j]
    b3.M2.out[,j] <- b3.M2.out[,1] + b3.M2.out[,j]
}
colnames(a3.M1.out)[1] <- colnames(a3.pairwise.out)[1]
colnames(a3.M2.out)[1] <- colnames(a3.pairwise.out)[1]
colnames(b3.M1.out)[1] <- colnames(a3.pairwise.out)[1]
colnames(b3.M2.out)[1] <- colnames(a3.pairwise.out)[1]


M <- nrow(a3.M1.out)

## permute the MCMC output
perm.inds <- sample(1:M, size=M, replace=FALSE)
a3.M1.out.perm <- a3.M1.out[perm.inds,]
perm.inds <- sample(1:M, size=M, replace=FALSE)
a3.M2.out.perm <- a3.M2.out[perm.inds,]
perm.inds <- sample(1:M, size=M, replace=FALSE)
b3.M1.out.perm <- b3.M1.out[perm.inds,]
perm.inds <- sample(1:M, size=M, replace=FALSE)
b3.M2.out.perm <- b3.M2.out[perm.inds,]
perm.inds <- sample(1:M, size=M, replace=FALSE)
a3.pairwise.out.perm <- a3.pairwise.out[perm.inds,]
perm.inds <- sample(1:M, size=M, replace=FALSE)
b3.pairwise.out.perm <- b3.pairwise.out[perm.inds,]


## storage objects
D.spear.a3.M1.holder <- rep(NA, M)
D.spear.a3.M2.holder <- rep(NA, M)
D.spear.b3.M1.holder <- rep(NA, M)
D.spear.b3.M2.holder <- rep(NA, M)
D.spear.a3.pair.holder <- rep(NA, M)
D.spear.b3.pair.holder <- rep(NA, M)

D.kend.a3.M1.holder <- rep(NA, M)
D.kend.a3.M2.holder <- rep(NA, M)
D.kend.b3.M1.holder <- rep(NA, M)
D.kend.b3.M2.holder <- rep(NA, M)
D.kend.a3.pair.holder <- rep(NA, M)
D.kend.b3.pair.holder <- rep(NA, M)


## M1 vs. pairwise
for (iter in 1:M){


    D.spear.a3.M1.holder[iter] <- Spearman(rank(a3.M1.out[iter,]),
                                           rank(a3.M1.out.perm[iter,]),
                                           48, 48)
    D.spear.a3.M2.holder[iter] <- Spearman(rank(a3.M2.out[iter,]),
                                           rank(a3.M2.out.perm[iter,]),
                                           48, 48)
    D.spear.b3.M1.holder[iter] <- Spearman(rank(b3.M1.out[iter,]),
                                           rank(b3.M1.out.perm[iter,]),
                                           48, 48)
    D.spear.b3.M2.holder[iter] <- Spearman(rank(b3.M2.out[iter,]),
                                           rank(b3.M2.out.perm[iter,]),
                                           48, 48)
    D.spear.a3.pair.holder[iter] <- Spearman(rank(a3.pairwise.out[iter,]),
                                             rank(a3.pairwise.out.perm[iter,]),
                                             48, 48)
    D.spear.b3.pair.holder[iter] <- Spearman(rank(b3.pairwise.out[iter,]),
                                             rank(b3.pairwise.out.perm[iter,]),
                                             48, 48)

    
    D.kend.a3.M1.holder[iter] <- Kendall2Lists(rank(a3.M1.out[iter,]),
                                               rank(a3.M1.out.perm[iter,]),
                                               48, 48, 48)
    D.kend.a3.M2.holder[iter] <- Kendall2Lists(rank(a3.M2.out[iter,]),
                                               rank(a3.M2.out.perm[iter,]),
                                               48, 48, 48)
    D.kend.b3.M1.holder[iter] <- Kendall2Lists(rank(b3.M1.out[iter,]),
                                               rank(b3.M1.out.perm[iter,]),
                                               48, 48, 48)
    D.kend.b3.M2.holder[iter] <- Kendall2Lists(rank(b3.M2.out[iter,]),
                                               rank(b3.M2.out.perm[iter,]),
                                               48, 48, 48)
    D.kend.a3.pair.holder[iter] <- Kendall2Lists(rank(a3.pairwise.out[iter,]),
                                                 rank(a3.pairwise.out.perm[iter,]),
                                                 48, 48, 48)
    D.kend.b3.pair.holder[iter] <- Kendall2Lists(rank(b3.pairwise.out[iter,]),
                                                 rank(b3.pairwise.out.perm[iter,]),
                                                 48, 48, 48)
    
}



cat("\n\n@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@\n")
cat("a3 - b3 precision (black vs. white coders)\n")
cat("-----------------------------------------------------------\n")
cat("                        D.spear                 D.kend\n")
cat("-----------------------------------------------------------\n")
cat("a3.M1                   ", round(mean(D.spear.a3.M1.holder), 3),
    "               ",
    round(mean(D.kend.a3.M1.holder), 3), "\n")
cat("b3.M1                   ", round(mean(D.spear.b3.M1.holder), 3),
    "               ",
    round(mean(D.kend.b3.M1.holder), 3), "\n")
cat("a3.M2                   ", round(mean(D.spear.a3.M2.holder), 3),
    "               ",
    round(mean(D.kend.a3.M2.holder), 3), "\n")
cat("b3.M2                   ", round(mean(D.spear.b3.M2.holder), 3),
    "               ",
    round(mean(D.kend.b3.M2.holder), 3), "\n")
cat("a3.pairwise             ", round(mean(D.spear.a3.pair.holder), 3),
    "               ",
    round(mean(D.kend.a3.pair.holder), 3), "\n")
cat("b3.pairwise             ", round(mean(D.spear.b3.pair.holder), 3),
    "               ",
    round(mean(D.kend.b3.pair.holder), 3), "\n")

cat("@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@\n\n")


D.hat.data <- data.frame(D.spear.a3.M1=D.spear.a3.M1.holder,
                         D.spear.a3.M2=D.spear.a3.M2.holder,
                         D.spear.b3.M1=D.spear.b3.M1.holder,
                         D.spear.b3.M2=D.spear.b3.M2.holder,
                         D.spear.a3.pair=D.spear.a3.pair.holder,
                         D.spear.b3.pair=D.spear.b3.pair.holder,
                         D.kend.a3.M1=D.kend.a3.M1.holder,
                         D.kend.a3.M2=D.kend.a3.M2.holder,
                         D.kend.b3.M1=D.kend.b3.M1.holder,
                         D.kend.b3.M2=D.kend.b3.M2.holder,
                         D.kend.a3.pair=D.kend.a3.pair.holder,
                         D.kend.b3.pair=D.kend.b3.pair.holder)


save(D.hat.data, file="MM3.2.2.D.hat.a3-b3.data.Rda")
                         
                         



########################################################################
## END a3-b3 comparisons
########################################################################






########################################################################
## START a4-b4 comparisons
########################################################################


## load MCMC output
load("MM3.2.a4.M1.out.Rda")
a4.M1.out <- a4.M1.out[sample(1:nrow(a4.M1.out), nrow(a4.M1.out),
                              replace=FALSE), ]
load("MM3.2.a4.M2.out.Rda")
a4.M2.out <- a4.M2.out[sample(1:nrow(a4.M2.out), nrow(a4.M2.out),
                              replace=FALSE), ]
load("MM3.2.a4.pairwise.out.Rda")
a4.pairwise.out <- a4.pairwise.out[sample(1:nrow(a4.pairwise.out),
                                          nrow(a4.pairwise.out),
                              replace=FALSE), ]

load("MM3.2.b4.M1.out.Rda")
b4.M1.out <- b4.M1.out[sample(1:nrow(b4.M1.out), nrow(b4.M1.out),
                              replace=FALSE), ]

load("MM3.2.b4.M2.out.Rda")
b4.M2.out <- b4.M2.out[sample(1:nrow(b4.M2.out), nrow(b4.M2.out),
                              replace=FALSE), ]
load("MM3.2.b4.pairwise.out.Rda")
b4.pairwise.out <- b4.pairwise.out[sample(1:nrow(b4.pairwise.out),
                                          nrow(b4.pairwise.out),
                              replace=FALSE), ]


## make MCMC output comparable
a4.M1.out <- a4.M1.out[,1:48]
a4.M2.out <- a4.M2.out[,1:48]
b4.M1.out <- b4.M1.out[,1:48]
b4.M2.out <- b4.M2.out[,1:48]

colnames(a4.M1.out) <- gsub("as.factor\\(photoID\\)", "", colnames(a4.M1.out))
colnames(a4.M2.out) <- gsub("as.factor\\(photoID\\)", "", colnames(a4.M2.out))
colnames(b4.M1.out) <- gsub("as.factor\\(photoID\\)", "", colnames(b4.M1.out))
colnames(b4.M2.out) <- gsub("as.factor\\(photoID\\)", "", colnames(b4.M2.out))
colnames(a4.pairwise.out) <- gsub("theta\\.", "", colnames(a4.pairwise.out))
colnames(b4.pairwise.out) <- gsub("theta\\.", "", colnames(b4.pairwise.out))

for (j in 2:48){
    a4.M1.out[,j] <- a4.M1.out[,1] + a4.M1.out[,j]
    a4.M2.out[,j] <- a4.M2.out[,1] + a4.M2.out[,j]
    b4.M1.out[,j] <- b4.M1.out[,1] + b4.M1.out[,j]
    b4.M2.out[,j] <- b4.M2.out[,1] + b4.M2.out[,j]
}
colnames(a4.M1.out)[1] <- colnames(a4.pairwise.out)[1]
colnames(a4.M2.out)[1] <- colnames(a4.pairwise.out)[1]
colnames(b4.M1.out)[1] <- colnames(a4.pairwise.out)[1]
colnames(b4.M2.out)[1] <- colnames(a4.pairwise.out)[1]


M <- nrow(a4.M1.out)

## permute the MCMC output
perm.inds <- sample(1:M, size=M, replace=FALSE)
a4.M1.out.perm <- a4.M1.out[perm.inds,]
perm.inds <- sample(1:M, size=M, replace=FALSE)
a4.M2.out.perm <- a4.M2.out[perm.inds,]
perm.inds <- sample(1:M, size=M, replace=FALSE)
b4.M1.out.perm <- b4.M1.out[perm.inds,]
perm.inds <- sample(1:M, size=M, replace=FALSE)
b4.M2.out.perm <- b4.M2.out[perm.inds,]
perm.inds <- sample(1:M, size=M, replace=FALSE)
a4.pairwise.out.perm <- a4.pairwise.out[perm.inds,]
perm.inds <- sample(1:M, size=M, replace=FALSE)
b4.pairwise.out.perm <- b4.pairwise.out[perm.inds,]


## storage objects
D.spear.a4.M1.holder <- rep(NA, M)
D.spear.a4.M2.holder <- rep(NA, M)
D.spear.b4.M1.holder <- rep(NA, M)
D.spear.b4.M2.holder <- rep(NA, M)
D.spear.a4.pair.holder <- rep(NA, M)
D.spear.b4.pair.holder <- rep(NA, M)

D.kend.a4.M1.holder <- rep(NA, M)
D.kend.a4.M2.holder <- rep(NA, M)
D.kend.b4.M1.holder <- rep(NA, M)
D.kend.b4.M2.holder <- rep(NA, M)
D.kend.a4.pair.holder <- rep(NA, M)
D.kend.b4.pair.holder <- rep(NA, M)


## M1 vs. pairwise
for (iter in 1:M){


    D.spear.a4.M1.holder[iter] <- Spearman(rank(a4.M1.out[iter,]),
                                           rank(a4.M1.out.perm[iter,]),
                                           48, 48)
    D.spear.a4.M2.holder[iter] <- Spearman(rank(a4.M2.out[iter,]),
                                           rank(a4.M2.out.perm[iter,]),
                                           48, 48)
    D.spear.b4.M1.holder[iter] <- Spearman(rank(b4.M1.out[iter,]),
                                           rank(b4.M1.out.perm[iter,]),
                                           48, 48)
    D.spear.b4.M2.holder[iter] <- Spearman(rank(b4.M2.out[iter,]),
                                           rank(b4.M2.out.perm[iter,]),
                                           48, 48)
    D.spear.a4.pair.holder[iter] <- Spearman(rank(a4.pairwise.out[iter,]),
                                             rank(a4.pairwise.out.perm[iter,]),
                                             48, 48)
    D.spear.b4.pair.holder[iter] <- Spearman(rank(b4.pairwise.out[iter,]),
                                             rank(b4.pairwise.out.perm[iter,]),
                                             48, 48)

    
    D.kend.a4.M1.holder[iter] <- Kendall2Lists(rank(a4.M1.out[iter,]),
                                               rank(a4.M1.out.perm[iter,]),
                                               48, 48, 48)
    D.kend.a4.M2.holder[iter] <- Kendall2Lists(rank(a4.M2.out[iter,]),
                                               rank(a4.M2.out.perm[iter,]),
                                               48, 48, 48)
    D.kend.b4.M1.holder[iter] <- Kendall2Lists(rank(b4.M1.out[iter,]),
                                               rank(b4.M1.out.perm[iter,]),
                                               48, 48, 48)
    D.kend.b4.M2.holder[iter] <- Kendall2Lists(rank(b4.M2.out[iter,]),
                                               rank(b4.M2.out.perm[iter,]),
                                               48, 48, 48)
    D.kend.a4.pair.holder[iter] <- Kendall2Lists(rank(a4.pairwise.out[iter,]),
                                                 rank(a4.pairwise.out.perm[iter,]),
                                                 48, 48, 48)
    D.kend.b4.pair.holder[iter] <- Kendall2Lists(rank(b4.pairwise.out[iter,]),
                                                 rank(b4.pairwise.out.perm[iter,]),
                                                 48, 48, 48)
    
}



cat("\n\n@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@\n")
cat("a4 - b4 precision (black vs. white distractor photos)\n")
cat("-----------------------------------------------------------\n")
cat("                        D.spear                 D.kend\n")
cat("-----------------------------------------------------------\n")
cat("a4.M1                   ", round(mean(D.spear.a4.M1.holder), 3),
    "               ",
    round(mean(D.kend.a4.M1.holder), 3), "\n")
cat("b4.M1                   ", round(mean(D.spear.b4.M1.holder), 3),
    "               ",
    round(mean(D.kend.b4.M1.holder), 3), "\n")
cat("a4.M2                   ", round(mean(D.spear.a4.M2.holder), 3),
    "               ",
    round(mean(D.kend.a4.M2.holder), 3), "\n")
cat("b4.M2                   ", round(mean(D.spear.b4.M2.holder), 3),
    "               ",
    round(mean(D.kend.b4.M2.holder), 3), "\n")
cat("a4.pairwise             ", round(mean(D.spear.a4.pair.holder), 3),
    "               ",
    round(mean(D.kend.a4.pair.holder), 3), "\n")
cat("b4.pairwise             ", round(mean(D.spear.b4.pair.holder), 3),
    "               ",
    round(mean(D.kend.b4.pair.holder), 3), "\n")

cat("@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@\n\n")


D.hat.data <- data.frame(D.spear.a4.M1=D.spear.a4.M1.holder,
                         D.spear.a4.M2=D.spear.a4.M2.holder,
                         D.spear.b4.M1=D.spear.b4.M1.holder,
                         D.spear.b4.M2=D.spear.b4.M2.holder,
                         D.spear.a4.pair=D.spear.a4.pair.holder,
                         D.spear.b4.pair=D.spear.b4.pair.holder,
                         D.kend.a4.M1=D.kend.a4.M1.holder,
                         D.kend.a4.M2=D.kend.a4.M2.holder,
                         D.kend.b4.M1=D.kend.b4.M1.holder,
                         D.kend.b4.M2=D.kend.b4.M2.holder,
                         D.kend.a4.pair=D.kend.a4.pair.holder,
                         D.kend.b4.pair=D.kend.b4.pair.holder)


save(D.hat.data, file="MM3.2.2.D.hat.a4-b4.data.Rda")
                         
                         



########################################################################
## END a4-b4 comparisons
########################################################################



sink()
